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Abstract. The statistics of multi- field inflation are investigated using the stochastic 
approach. We analytically obtain the probability distribution function of fields with 
| the scaling approximation by extending the previous work by Amendola. The non- 

Gaussian nature of the probability distribution function is investigated decomposing 
the fields into the adiabatic and isocurvature components. We find that the non- 
pT) ' Gaussianity of the isocurvature component can be large compared with that of the 

■ adiabatic component. The adiabatic and isocurvature components may be correlated 

'. at nonlinear order in the skewness and kurtosis even if uncorrelated at linear level. 
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1. Introduction 

The inflationary universe is now the standard paradigm that explains the observed 
universe naturally QjJ One of the most important predictions of inflation is the 
primordial fluctuations from which the large scale structure of the universe originates. 
Many models of the inflation have been proposed. Multi-field inflation is a natural 
alternative model to the standard single field inflaton model (e.g., EJ El El E] ) • In 
such the multi-field inflation model, the isocurvature perturbation arises in addition to 
the adiabatic one, though the evolution of the isocurvature perturbations is model- 
dependent. Testing predictions of these inflation models will be the subject of 
the (future) precision cosmology by confronting with observation. Statistics of the 
primordial fluctuations is very important in such the test of the inflation models. In 
particular, the non-Gaussian nature of the primordial fluctuations is considered to be 
an important clue [S], which stimulates the recent investigations on the higher order 
cosmological perturbations [IH EH HH H21 EH CH HS1 ESI HZl HH1 HSl EH] - 

The stochastic inflation was developed by Starobinsky and many works have been 
done |2U E21 |2H1 Ell- The stochastic approach to inflation is useful in investigating 
the statistics of the fluctuations [2H1I2E]- In the stochastic inflation, the field dynamics 
coarse-grained over the horizon size is described by the Langevin equation, and the 
field value of the inflaton is the probability variable. Instead of solving the Langevin 
equation, we may solve the corresponding Fokker-Planck equation, which describes the 
evolution of the probability function P{4>)- It is known that the scaling approximation 
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is useful to describe the evolution of the probability function during the slow-roll period 
[2*7] . With the use of the scaling approximation, the probability function can be obtained 
analytically [2*""j . The stochastic approach to the multi- field inflation has been considered 
in j2I3 EOl UH]- The probability function for the multi- field inflation was constructed 
using the scaling approximation in an analytic manner by Amendola [SI]. However, 
the constructed probability function does not have the symmetry in exchanging the 
field parameters. Such the symmetry is desired for a consistency of the theory. As an 
extension of this previous work, we develop a useful analytic formula for the probability 
function of the multi-field stochastic inflation using the scaling approximation. 

The organization of the present paper is as follows: In section 2, we summarize 
the stochastic approach to the multi-field inflation. Then we obtain the probability 
distribution function using the scaling approximation as an extension of the work by 
Amendola ?H"|. In section 3, based on the probability function, we investigate the non- 
Gaussian feature of the scalar fields during the inflation by decomposing the scalar 
fields into the adiabatic and the isocurvature components. In section 4, we mention 
about connecting our result with observation. Section 5 is devoted to summary and 
conclusions. Throughout this paper we use the unit c = h = 1, and follow the convention 
(+,-,-,-)■ 



2. Formulation 



We consider the chaotic inflation model with the two scalar fields and x with the 
Lagrangian density 

c = ivvv M - Vttf) + -v^v^ - V 2 (x), (1) 

where V\{4>) and V^x) are the potential, 

W) = ^ 4 (rf • (2) 




ff 2 = " w ~7 2 "' , (4) 



^(x) = 7r~ a ( _ ) • (3) 



where Ai and A2 are the non-dimensional coupling constants, n is the parameter, and 
we introduced a~ 2 = 8ttG. During the slow-role regime, neglecting the kinetic terms, 
the Friedman equation gives 

v 1 (<f>)_+y 2 (x) 

3a 2 

where H = a(t)/a(t) is the Hubble parameters, a(t) is the scale factor, the dot denotes 
the derivative with respect to the cosmic time t. 

We adopt the stochastic approach to the inflation. The basic equation, the Langevin 
equation, is obtained by coarse-graining the field, averaging over the short wave length 
modes of the evolution equation of the field operator. For the two-field inflation, we 
have 

1 -Vi„ 



3H +m ^ 
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X = -jjf- + 9Vx: (6) 

where 

H 3/2 

9 = -^p, (7) 

and rj^ and r] x describes the stochastic noises, which have the nature of the Morkovian 
Gaussian noise with zero mean and the correlation 

(vM) = 0, (8) 
(rH(t)v j ^)) = S iS ^S{t-f), (9) 

where r/j and rjj denote rj^ and rj x . Thus the long-range scalar fields during the slow-roll 
regime are described by the Langevin equations with the classical drift force — V\, 4,/SH 
or — V2, and the noise whose amplitude is in proportion to g. Note that, even if the 
field potentials consist of independent function of and Xi ^ ne dynamics is coupled to 
each other because the Hubble parameter H depends on the both fields. Note that 
we have assumed (T)<p{t)T) x (t')) = 0, based on the assumption that the fine-grained 
component of the field <fi and x can De treated as free in the leading approximation 
|29| I3T] . Then, we may assume the Langevin equations for these field with the Gaussian 
noise with the correlation (^(t)^(i')) = ^ — 0> i n a similar way to 

the single field stochastic inflation. However, this assumption might be oversimplified 
for the models with general coupling between and %. 

Here let us consider the background solution for the Langevin equation without the 
noise term, 



3H 



(10) 



From these zero-noise field equations, we have 



(12) 



dxd ^xT 11 

which gives the attractor solution, assuming that the classical paths converge to it after 
a short transient, 

d = L^-^Xd, (13) 

where L = \ 2 /\\. In the present paper, we refer to the attractor solution as the 
classical solution. We will solve the Fokker-Planck equation corresponding to the 
Langevin equation, along the classical attractor trajectory. For the latter convenience, 
we introduce the parameter 9 by 

cos6= h = (1 + £l/(l-n))-l/2 (14) 
sin = X* = ( L V(n-l) + !)-l/2_ (15) 
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Our formulation relies on rewriting the Langevin equation in term of the variables 
r and G, introduced in the below, instead of <fi and \- First, we define the variable r by 

r = <p cos 6 + x sm &■ (16) 
From the Langevin equation for and \i we niay write 
—Act 3 /r\ 2 " -1 

r = ^-(-J + #(cos%, + sin#r/ x ), (17) 
where H should be regarded as 

1 \rr 2 / r \ 2n 

H 2 = —r (Vi(rcos0) + V 2 (rsm6)) = ^- ( - ) (18) 

with defined 

A = Ai cos 2 " 6 + A 2 sin 2 " 6 = (A} /(1 ~ n) + A^ 1 ^) 1 -" (19) 

The approximation of equation (fTSj) is similar to that in [3T]. In the strict sense, H 
is the function of and x> t> u ^ m equation (fTHj) . if is approximated as H((f>,x) = 
H(4> c i(r),Xci(r)). The other variable 6 is defined as 

^(l-n) ^ X 2(l-n) 



9 2(n-l) 
It is easy to show that satisfies 




A 2 \<J 



(20) 



e= ^-^Ucr 2 . (2i) 

We consider (JT7J) and (|2*T|) . instead of the original Langevin equations (0) and (JHJ). 

It will be important to mention about the physical meanings of the variable r and 0. 
It is clear that r is the variable parallel to the classical trajectory. On the other hand O 
can be regarded to represent the component perpendicular to the classical trajectory. It 
is worthy to note that the cross correlation of the noise terms in the Langevin equations 
r and O is zero, i.e., 

^ (cos 61^ + sin%) - ^) ) = °" ( 22 ) 

Now we consider solving the Fokker-Planck equation corresponding to the Langevin 
equations for r (and O). To find the solution analytically, we adopt the scaling 
approximation, which was done in reference [2H] for single inflaton models. In the 
slow-roll regime of inflation, the drift force dominates the dynamical evolution of the 
scalar fields, and the quantum noise is small. Then the coarse-grained probability 
distribution develops a sharp peak around the classical slow-roll solution. This is the 
basic background that the scaling approximation may be used to describe the behavior 
of the probability distribution in the slow-roll regime. The Langevin equation for r is 
similar to that for the single field inflaton model, and we follow the reference j2H] to find 
the analytic solution with the technique of the scaling approximation. In Appendix A, 
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the derivation of the probability function is described. The analytic probability function 
can be written as function of (f) and \ as 



Psc{<f>,X) 



exp 



2v u 



where J(<p, x) is the Jacobian of the transformation 
and the variances are obtained as 

i nr 



2v 



22 



\J{<f>,X)\, (23) 



(24) 



and 



I67H V 2nA V o 



A 



1- — 



r 



fn > 2) 



(25) 



V22 



+ 



192tt 2 I A 2 cos 6 fl A|sin 6 ^ 



-In — 



r 



1 



967r 2 n 2 (n - 2) 

' y .\ 2(2— n) 



A? cos^ 2 ™- 1 ^ 

' r \ 2(2-n)" 
(7 



A 2 sin 2(2n-1)^ 



(n 



fn > 2) 



(26) 



and £ is related to r(0, x) by equation (jA.6|) and £* = £(r c i). 

Here we mention about the boundary condition of the probability distribution 
function. The reflecting boundary condition, with which the normalization of the 
probability is preserved, is adopted in usual. In this case, other reflective term should 
be added to the probability function. However, as long as we are concerning with the 
behavior of the probability function around the peak of the classical slow-roll solution, 
we may work with the solution ()23j) practically. 

Figure 1 shows the contours of the probability function to make an example of 
its time evolution. At the earlier time (left panel), the probability function is almost 
spherically distributed. At the middle time (center panel), the probability function 
becomes rather broad, then it becomes elongated along the direction of the classical 
path at the latter time (right panel). The behavior depends on the model parameters 
as well as the initial condition. 



3. Non-Gaussian behavior of the probability function 

In this section we investigate the evolution of the probability function (}23|) . focusing on 
the non-Gaussian statistical nature. We start from introducing the variable Ar and As 
by 

Ar \ = / cos 6 sm.6 \ f <j) - 4> d \ 
As J { - sin9 cos9 J \ X - Xd J ' 

We may regard Ar (As) as the component that contributes to the adiabatic 
(isocurvature) perturbation IE] • We rewrite the probability distribution function in 
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terms of Ar and As. The expansion can be done in a straightforward manner. But the 
expression is rather long, then the result is given in Appendix B. 

Here let us introduce the e-folding number from a time t to the final time of the 
inflation tf, 



N(t) = 

where tf is defined by Vi(4> c i 



H(t')dt', 

' v 2 ( Xd ) 



(28) 



^ + Xli) A which gives r cl (t f ) 2 = 2nV 2 /3. 
Instead of using t (or r c i(t)), we use N as the time variable. The classical solution is 
related to N by 

2 



rci{t) 



n + 6N 



J\f < 1 



(29) 



\ r* J n + 6iV* 

for n > 2, where N* = N(t*). 

With the use of Af (or N), the part of the probability function of £, equation (jB.lj) 
is written 



exp 



2vn 



cxp 



71-1 + 



n + 2 



V r d 



+ 



'(n-l)(7n-3) , 3(n - l)(n +_2)^_„ 



12 



+ 



4 



+ (n + 2)(30-n) A , 2(2 _ n) ^| fAr 



48 



and R„ is 



96tt 2 



(|-) n n 2 (Ar 2 - !)-!(„ + 6iV)~ n 



(30) 



(31) 



from (jB.2|) . On the other hand, the part of the probability function of is also written 
in the similar way, i.e., equation (jB.3|) . with 

9 



_ 96tt 2 

On — ~ X 

A 



O)2 ;i + 3A0- 2 (-ln./vT 1 
(n - 2)n 2 / 3 ^ n 
2 V2^ 



(71 



n + 6iV)~ n (l - A/" ) (n>2) 



From equations (j3*T|) and 
S " _ J 2 



we easily have 
(A"~ 2 — 1)(— lnATr 1 



(n = 2) 



71 



■(JV- 2 -!)^-^ 2 )- 1 (ti > 2) 



(32) 



(33) 



and S n /R n > 1 (see Figure 2). Because the variance of Ar and As is in proportion to 
l/R n and 1/S n , respectively, then it can be read that the amplitude of the fluctuation 
of the isocurvature component is smaller than that of the adiabatic component. This is 
consistent with the behavior of the probability function shown in Figure 1. 
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We next consider the deviation of the probability function from the Gaussian 
statistics, which is expressed by the terms in proportion of Ar 3 , As 3 , and the other 
higher order terms. To investigate the deviation of the probability distribution function 
from the Gaussian statistics quantitatively, we first define the ratio 
P, c ( + A r , + A 3 ) 

P, c (-Ar,-As) V ' 

Thus T represents the asymmetry of the probability function, which can be related to 
the dimensionless skewness (see below). In the limit of the Gaussian statistics, T — 1 
approaches to zero. This is a very simple but useful estimator of the local skewness 
around the Gaussian peak introduced in reference [321 • From our probability function, 
we have the expression of the leading order 

T - 1 »2\(n - 1 + ^A/-< 2 -" 2 ) R» (f-^ - (2n - l)r„S„ (^ 

+ 2 ( 2n - 1)5 „(^)(^)], (35, 

for n > 2, where r„ is defined by (jB.5|) . Our result of the adiabatic component with 
n = 2 reproduces the result in reference [3*2] . 

We can read a few characteristic features of the skewness from expression (|33|) . 
First, the evolution of the skewness, in general, depends on the parameters of the 
potential. It also depends on the initial time N(t*), which we have adopted to solve 
the Fokker-Planck equation. This is unavoidable in our scheme. Let us consider the 
skewness of the adiabatic component, the first term of the right hand side (r.h.s.) of 
()35|). To estimate the skewness quantitatively, we consider the first term at the one 
sigma deviation with respect to Ar, with replacing Ar/r d by 1/V2i? n , 

„ a= („_ 1 + (36) 

In the similar way, considering the second term of (r.h.s.) of ()35|) at the one sigma 
deviation with respect to As, with replacing As/r c / by l/\ / 2S n , we define 

A 5 = -(2n-l)r n -^=. (37) 

Figure 3 plots Ar and As as function of iV for the models Ai = A2 = 1CT 12 , and 
Ai = 1CT 12 and A2 = 1CT 13 . For the model Ai = A2, we have As = because T n = 0. 
Here we extensively used our result to the case n = 1.5. We see that As is larger than 
Ar for the model Ai = 10~ 12 and A2 = 10~ 13 . Figure 4 plots the same as Figure 3 
but with n = 2. Similarly, Figures 5 and 6 are same as Figure 3 but with n = 3 and 
n = 4, respectively. From Figures 3-6, we see that the non-Gaussian estimator Ar and 
As becomes larger as n becomes larger, in general, depending on A. But the ratio of As 
to Ar becomes large as n becomes small. This comes from the fact that the skewness of 
the isocurvature component As is in proportion to Y n . Figure 7 plots the ratio As/Ar 
at a fixed time as function of L for various models of n. This result demonstrates that 
the skewness of the isocurvature component is independent of the adiabatic component, 
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and can be large compared to the adiabatic one. The asymmetry in the two potentials 
enhances the non-Gaussianity of the isocurvature component. In reference j^Sl, the 
importance of the non-Gaussianity of the isocurvature component was discussed in a 
different context of the investigation of multi-field inflation. 

The third term of the r.h.s. of represents the mixing of the isocurvature and the 
adiabatic components. This suggests that the adiabatic and isocurvature components 
may have a correlation at nonlinear order in general, even if it have no correlation at 
the linear order [34*] . 

Similar to the skewness , the terms in proportion to Ar 4 and As 4 represent the local 
kurtosis of the probability function. It can be read that the coefficient of As 4 depends 
on T n but the coefficient of Ar 4 does not. This means that, similar to the case of the 
skewness, the kurtosis of the isocurvature component is independent of the adiabatic 
component and can be large compared to it. 



4. Connecting the result with observation 



Here let us briefly mention about the detectability of the skewness and kurtosis. As 
mentioned in section 1, many works related to the non-Gaussian nature of the primordial 
fluctuations have been done. The useful quantity f^L is introduced to characterize the 
non-Gaussianity of the gravitational potential $, by expanding it as 35J 

$ = <P + f N L& 2 - (f 2 )) + OUnl ), (38) 

where <p is a zero-mean Gaussian probability variable. Here assume <p and $ are the 
variables with one degree of freedom, it is easy to construct the probability distribution 
function of $, 

$ 2 - 2f NL &~ 



P($) oc exp 



2 ($ 2 ) 



(39) 



neglecting the terms proportional to f^i ■ X Then we define the local skewness for the 
probability function P($), in the similar way to the previous section, by 

Pi 



<$ 2 » 



<$ 2 » 



-1-2/ 



NL 



$ 2 ) 1/2 . 



(40) 



We can write 7$ — 1 in terms of the dimensionless skewness defined by ($ 3 ) / ($ 2 ) 3 ^ 2 . 
Approximating the probability function as 



we have 



exp 



($3) 



$ 2 - 2f NL & 
2($ 2 ) 



1 + /. 



$ 3 



NL 



<$ 2 >, 



exp 



$5 



2($ 2 ) 



15W($ 2 > 



15 



(% - 1). 



(41) 



(42) 



($2)3/2 -—V v- / 2 

| In deriving this expression we redefined + /nl (^& 2 ) by $ to eliminate the term proportional to $ 
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In general, the evolution of fluctuations depends on each model. However, we 
here introduce the assumption that the dimensionless skewness preserves, that is 
|7$ — 1 1 ~ \T — 1|, even if the amplitude of the fluctuations may change during the 
evolution of the universe. Based on this assumption our result can be connected with 
observation. In this case the predicted local skewness is related to f^L as 

\T-l\ = 2f NL J@*) or lLz3 = 2f NL . (43) 

v<$ a > 

If we adopt J ($ 2 ) = 10~ 5 , the vertical axis of Figures 3-6 can be read as IJnl- 

The detectable minimum value of Jnl is discussed in references [33], f^L = 20 and 
5 for WMAP and Planck data, respectively. Thus the non-Gaussianity in our model 
with n > 3 might be detectable. However, this is based on the rough estimation because 
we have not strictly considered the observational constraint on the model parameters 
from the spectrum index and the amplitude of the perturbation. For definite conclusions 
we need incorporate these constraints. 

Similar to the skewness, we can estimate the local kurtosis, being of order C(($ 2 )). 
This is smaller than the skewness by a factor of 10 -5 . 



5. Summary and conclusions 

In the present paper, we have investigated the statistical nature of the multi-field 
inflation with the stochastic approach, in which the evolution of the scalar fields 
is described by the Langevin equations. To solve the corresponding Fokker-Planck 
equation analytically, we adopt the technique of the scaling approximation. This was 
first considered in We have improved the previous investigation requiring the 

probability function to have the symmetry in exchanging the field parameters. We have 
investigated the non- Gaussian nature of the probability function by decomposing the 
fields into the adiabatic and the isocurvature components. 

We have introduced the local skewness as an estimator of the non-Gaussianity, and 
have found the followings: (1) The amplitude of the skewness is determined by the 
combination of the coupling constant Ai, A 2 and n. The amplitude of the skewness of 
the adiabatic component is determined by the combination of the coupling constant A 
and n, while that of the isocurvature component depends on A, n and the ratio of the 
coupling constant T n . Then the skewness of the isocurvature component can be large 
compared with that of the adiabatic one due to the factor T n . (2) The similar feature 
can be seen in the kurtosis of the adiabatic and isocurvature components. (3) The 
adiabatic and isocurvature components may have correlation at the nonlinear order in 
the skewness and kurtosis, even if there is no correlation at the linear order. 

We have considered the connection between our result and observation, based on the 
assumption that the dimensionless skewness of fluctuation is preserved. This argument 
shows that the skewness will be useful to constrain the potential parameter. However 
the assumption is not trivial and there remains room for more discussions for precise 
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comparison with observation. Furthermore, in the present paper, we have only worked 
with the simple chaotic inflation models with the power-low potentials. Investigation of 
multi-field inflation with more general potentials will be an interesting subject, including 
the confrontation with cosmic microwave anisotropy observation. Testing the robustness 
of the approximation adopted here will also be required as a future investigation. 
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Appendix A. Scaling transformation and probability function 



To find the scaling solution of equation (j!7j) . two transformation on the variable 
is performed to yield a Langevin equation with drift term vanished. The first 
transformation is defined by 

g[r') 

where is an arbitrary constant, but we here chose = oo to get z with the same 
sign of r(> 0). In terms of this new variable, the Langevin equation is 



la 7 ' 2 



r[z) 



2n-l 



(cos 9rj$ + sin 9r] x ). 



(A.2) 



Thus, with the first transformation, the multiplicative (r-dependent noise) process 
reduces to the additive (constant noise) process, and the variable z is referred to as 
the constant-diffusion variable. Then we perform the second mapping, requiring that 
the resultant Langevin equation has no drift term, 

z->£ = F- 1 [e At F(z)]. (A.3) 

The function F is found in the form, 

A 



F(z) = exp 



-dz 



where D z is the drift term of equation (jA.2|) . 



D, 



la 7 ' 2 



3H 5 / 2 (r(z)) 



r[z) 



a 



2n-l 



(A.4) 



(A.5) 



and A = — D z (z*) is the separation constant and z* is chosen so that dF(z*)/dz = 1. 
Through these transformation, we obtain 

1 /12\ 3 / 4 fr 



£(r) 



2 V A 
2 

3n - 2 

' r 



1 'rd\ 2 



x 



a 



A 

2-n 



a J 



2-n 



2-n 



(3n-2)/(2n-4) 



(n = 2) 



(n > 2) 



(A.6) 



where is the initial value of the classical solution, and we used the classical solution 
r c i(t), defined as 



r_d 



exp 



'Ml -at 



[n 



1 + (n - 2) 



(7 



n-2 



'2nA 



-at 



l/(2-n) 



(A.7) 



fn > 2) 



which is given by solving equation ()17j) neglecting the noise term. Now the Langevin 
equation for £ is 

(n+2)/(3n-2) 



(cos + sin 8r]^ 



(Ai 
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We follow the standard choice of the initial condition for the probability distribution, 
the delta-function distribution at the initial time. Furthermore, assuming the limit of 
small diffusion, we adopt the Gaussian form of the probability distribution function in 
the scaling approximation 



P(£) oc exp 



2v n 



(A.9) 



where v\\ is given by expression (|2o|) . and = £(r c z). 

On the other hand the Langevin equation for the variable G has no drift term 
initially, then we write the probability distribution function in the Gaussian form 
assuming the limit of small diffusion 

G 2 



P(G) oc exp 
where v 2 2 is given by (pjj) . 



2v 



22 



(A.10) 



Appendix B. Expansion of the probability function 



In this Appendix, we present the expression of the probability function expanded in 
terms of Ar and As. The variable £ is written in terms of r, equation (jA.6|) . Then the 
part of the probability function written by £ provides us with information about the 
adiabatic component, which can be expanded as follows, by setting r = r c \ + Ar, 

(£-£*) 2 " 



exp 



2v u 




exp 



f{n- l)(7n-3) 3(n - l)(ra + 2) /r d \ 2 ~ n 
+ V 12 + 4 \~J 

(n + 2)(30-n) /r c/ \ 2 ( 2 - n )\ /Ar 
+ 48 In J J W 



(B.1) 



for n > 2, where we defined 
96tt 2 , 



— — n 
X 



1 /^y(2-n) 



a J \ a ) V cr J 
On the other hand, the part of the probability function of G is expanded as 

02 



(B.2) 



exp 



2iJ 2 2 



exp 




As\ 2 . x fAs^ 3 

— + 2n-lT„ — 
r c i J V r d 



- 2(2n - 1 
2n 



Ar 



As; 2 



+ 



_L (( i4„-3)r.' + 8o )(0 



- (2n- l)(4n-l)T r 



Ar\ /As 



+ (2ra- l)(4n- 1) 



Ar\ 2 /As 
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(B-3) 



where we defined 



96tt 2 
&, = x < 
A 



-In 

n — 2)n 



2 



I- 1 /^N^-n) 
(J 



r \2(2-n) / r \2(2-n) 



(7 



1 -1 



a 



(n = 2) 

4(l-n) 



(B.4) 



(n > 2) 



and 



r n = tan0 - cot = (l/L) 1 ^"- 1 ) - L V2(n-i)_ 



(B.5) 
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Figure 1. Contours of the probability distribution function. Here the potential with 
Ai = 10~ 5 , A 2 = 10~ 4 , and n — 2 is considered. The initial time is chosen iV* = 100 
and the left, middle and right panels show the contour at the time of the c- folding 
N = 95, 60, and 20, respectively. (See equation (28) for the definition of N) 




Figure 2. The ratio S n /R n is plotted as function of N. The initial time is chosen 
-/V* = 100. The parameter of the potential n is chosen n = 1.1, 2, 3, 4 from bottom 
to top. Here we extensively used our result to the case n < 2. 
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Figure 4. Same as figure 3, but with n = 2. 
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Figure 7. The ratio As/Ar at the time N — 60 as function of L for the models 
n = 2, 3, 4 (from top to bottom). Here we fixed AiA 2 = 1CP 24 , and N* = 100 is 
adopted. 



